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As part of the UKQCD and RBC collaborations' TV/ = 2+1 domain- wall fermion phenomenology 
programme, we calculate the first two moments of the light-cone distribution amplitudes of the 
pseudoscalar mesons n and K and the (longitudinally-polarised) vector mesons p, K* and (j>. We 
obtain the desired quantities with good precision and are able to discern the expected quark-mass 
dependence of SU(3)-flavour breaking effects. An important ingredient of the calculation is the 
nonperturbative renormalisation of lattice operators using the RI'/MOM technique. 



I. INTRODUCTION 

Light-cone distribution amplitudes (DAs) are impor- 
tant nonperturbative quantities which (within the frame- 
work of collinear factorisation) parameterise in partonic 
terms the components of the hadronic wavefunction that 
control hard exclusive processes. Such processes provide 
hadron structure information complementary to that ob- 
tained from hard inclusive reactions. 

The structure functions for inclusive processes are 
more accessible both experimentally and theoretically 
owing to their larger cross sections and branching ratios, 
simpler final-state detection and more straightforward 
factorisation properties. They do not specify the phases 
and correlations which would constitute amplitude-level 
hadron structure information, but probe instead the 
bound states' partonic content. Deep-inelastic scattering 
processes, for example, are controlled by the charge and 
momentum of the struck parton and are insensitive to its 
relation to the other hadronic constituents. The associ- 
ated parton distribution functions (PDFs) are therefore 
single-particle probabilities, which reveal nothing about 
the role of particular Fock states or of correlations be- 
tween quarks and gluons. 

Distribution amplitudes, involved in exclusive pro- 
cesses, always appear in convolutions and, unlike the 
PDFs, are not directly measurable. These exclusive pro- 
cesses are dominated by specific partonic configurations. 
The outgoing quarks and gluons are unlikely to form a 
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given final-state hadron unless either they are approxi- 
mately collinear with small transverse separation, or one 
of the partons carries almost all of the hadron's momen- 
tum (the soft overlap or Feynman mechanism). In the 
former case, the basis for collinear factorisation hard 
gluon exchange must occur to allow the struck or decay- 
ing parton to communicate with the others, turning them 
to the final direction. Since more partons require more 
hard gluons, exclusive cross-sections and decay rates are 
dominated by the valence Fock state at leading-order in 
Q 2 , up to soft effects. 

Hard exclusive processes are therefore controlled at 
leading-order by the distribution amplitudes of leading- 
twist (an operator's twist is the difference between its 
dimension and its spin): essentially the overlap of the 
hadronic state with the valence Fock state in which, for 
a meson, the collinear quark-antiquark pair have small 
transverse separation and carry longitudinal momentum 
fractions u and u = 1 — u. The pion's electromagnetic 
form factor at large Q 2 , for example, can be written as 
a convolution of distribution amplitudes (f> v {u, Q 2 ) for 
the incoming and outgoing pions with a perturbatively- 
calculable hard-scattering kernel. Higher-twist DAs asso- 
ciated with power-suppressed contributions originate in, 
for example, higher Fock states [2|. We consider only 
leading, twist-2, DAs in this paper. 

The phenomenological importance of hard exclusive 
processes has grown since collinear factorisation was first 
established for cases such as the pion's electromagnetic 
form factor and the 77*71 transition form factor around 
30 years ago 0, 0-11] ■ Of particular note is the theoret- 
ical description of hadronic B decays, which have been 
studied in detail by BaBar and Belle and will be studied 
by LHCb and at super- B factories in order to constrain 
the CKM matrix and to understand CP violation. Fac- 
torisation is more difficult to establish in i?-physics be- 
cause the hard collinear and soft mechanisms contribute 
at the same order in I/77J&. Two approaches have been 
developed. In the QCD factorisation framework it has 
been shown that collinear factorisation can be applied to 
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leading order in l/m& to a large class of nonleptonic B- 
decays [HI]. Soft-collinear effective theory (SCET) 0- 
HTj | aims to provide a unified theoretical framework for 
the factorisation of both hard-collinear and soft effects. 
In both cases, distribution amplitudes play an important 
role as nonperturbative inputs in flavour physics. 

In this paper we focus on the distribution amplitudes of 
the light pseudoscalar and longitudinally-polarised vector 
mesons, since, as we shall discuss in Sec. II Bl their lowest 
moments are of phenomenological interest and are calcu- 
lable on the lattice. For pseudoscalars, these quantities 
are relevant for decays such as B — > irir and B — > kK; 
they also appear in light-cone sum rule (LCSR) expres- 
sions for the form factors of semileptonic decays such 
as B — > ixlv. For hard exclusive processes involving the 
light vector mesons p, K* and 0, polarisation-dependence 
can reveal much about the underlying dynamics, with 
the longitudinally- and transversely-polarised final vec- 
tor meson states often involving different aspects of weak 



12j. Examples are the exclusive 
rare radiative B — > p7 or non- 



interaction physics 
semileptonic B — > plvi 
leptonic, e.g. B — > np, decays of B-mesons, which are 
important for extracting CKM matrix elements. 



is the path-ordered exponential needed to maintain gauge 
invariance, p is a renormalisation scale, u is the momen- 
tum fraction of a quark, u — 1 — u and is the polar- 
isation vector for a vector meson with polarisation state 
A. The distribution amplitudes are normalised by 



du 4>(u, p) = 1 . 



(4) 



The definitions above involve the pion and rho-meson 
decay constants defined by 



|?27/i75Q'lK0)) = ifrrPp, 

|¥ 2 7a<9i|p(p; A)) = fptUpE^ 



(5) 
(G) 



The vector meson decay constant, / p , and its coupling to 
the tensor current, fj, are of interest in their own right 

and we have previously calculated [l3[ the ratios fy / fv, 
for V £ {/?, K* 7 4>} as part of our domain-wall fermion 
(DWF) phenomenology programme. 



Moments 



A. Definitions 

Mesonic light-cone DAs are defined from meson-to- 
vacuum matrix elements of quark- ant iquark light-cone 
operators, which are non-local generalisations of those 
used to define the decay constants. For example, for pi- 
ons 

(°|92( 2 )7^75^(^, -z)qt(-z)W(p))z*=o = 

iUVv I due^-^- z ^(u,p) (1) 
Jo 

and for longitudinally-polarised rho-mesons 
( I<72( Z )7^( Z , -z)qi{-z)\p{p; A)) z 2 =0 = 

f P m p p v / due^-^^p), (2) 

P ' z Jo 

where 

T(z, -z) = V cxp(-ig J dw^A^wfj (3) 
Zl I 



Moments of light-cone DAs are defined by: 

(f%0«)= / due 4>{u, p), (7) 
Jo 

where £ = u — u = 2u — 1 is the difference between the 
longitudinal momentum fractions. 

Since the moments arc obtained from matrix elements 
of local operators [14| we can study them using lattice 
QCD. The light-cone matrix elements which define the 
DAs themselves are not amenable to standard lattice 
techniques, since in Euclidean space the light-cone has 
been rotated to the complex direction. By expanding 
the non-local operators on the light cone, we obtain sym- 
metric, traceless twist-2 operators. With the following 
conventions for continuum covariant derivatives, 

— > — > <— <— <->•<——>• 

Dp = dp + igA p , Dp = dp- igA p , Dp = Dp- Dp, 

(8) 

the expressions relating the moments of DAs to the cor- 
responding local matrix elements are: 



(0\q(Qh P l5Dps(0)\K(p)) = (e)KfKPpPp , (9a) 

(0\q(0)l P -f5DpDM0)\Ap)) = -i(i 2 ),hP P PpP, , (9b) 

(0\q(0hpDp S (0)\K*( P ,X)) = (e) ll K *fK*m K ^ (p^W +p«4 A) ) , (9c) 

(0|«(0)7A^«(0)|p(p, A)) = f P m p ~ (s p X) Pp P „ + e^p v p p + e^ PpP p) . (9d) 
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The operators in the matrix elements above are all to be 
considered symmetric and traceless in the free Lorentz 
indices. Meson-meson rather than meson- vacuum ma- 
trix elements of the same operators lead to moments of 
generalised parton distributions (GPDs). 

Recent analyses, especially those based on QCD sum 
rules, deal instead with the Gegenbauer moments, which 
arise from a conformal expansion (l5l [l6j , in which the 
conformal invariance of (classical) massless QCD is used 
to separate longitudinal and transverse degrees of free- 
dom, analogous to the partial wave expansion in ordi- 
nary quantum mechanics. All dependence on the longi- 
tudinal momentum fractions is described by orthogonal 
polynomials that form an irreducible representation of 
the collinear subgroup of the conformal group, SL(2,M). 
The transverse-momentum dependence is represented as 
the scale-dependence of the relevant operators and is gov- 
erned by renormalisation-group equations. The different 
'partial waves', labelled by different conformal spins, do 
mix but not to leading-logarithmic accuracy. Conformal 
spin is thus a good quantum number in hard processes 
up to small corrections of order a 2 . 

The asymptotic Q 2 — > oo DA is known from pertur- 
bative QCD: </> as = 6uu. For the leading-twist quark- 
antiquark DAs that we are interested in, the conformal 
expansion can then be conveniently written as: 



oo 

u,fi) =6uu( 1 + Va„(/i)Cf (2u- 1) 



(10) 



3/2 

where C„ are Gegenbauer polynomials. To one-loop 
order the Gegenbauer moments, a n , renormalise multi- 
plicatively [l6[: 



are combinations of ordinary moments of equal and lower 
order, e.g.: 



a-2 



7 

12 



H 2 



i 



(13) 



C. Status 

In this section, we summarise what is currently known 
about leading-twist light-meson distribution amplitudes. 
For mesons of definite G-parity, there is a symmetry un- 
der the interchange « f> « of the two momentum frac- 
tions. In these cases, the distribution amplitude is an 
even function of £ = u — u and the odd moments there- 
fore vanish. Thus, (C 1 )^, and (C 1 )!} all vanish, while 

{^ 1 )k and (£ 1 )1 f . are SU(3)-flavour breaking effects. 

Since {S})k is the average difference between the frac- 
tions of longitudinal momentum carried by the strange 
and light quarks, 

(£V(M) = / du(2u-l)0 Jf (u > At) = <2u-l), (14) 
Jo 

we may expect from the constituent quark model that 
the sign of (^k — faf" is positive and this is indeed the 
case. In fact, is an important SU(3)-breaking pa- 

rameter and is relevant for predictions of B-decay transi- 
tions such as B — > K , K* 12311 . For example, a light-cone 
sum rule analysis leads to [24j : 



f± _ (0) fK n . K\ . 



(15) 



a n {n) = a n (no) 



(7(n)-7(o))/A) 



,a s (po)J 

The one-loop anomalous dimensions are: 



(11) 



7(n) = 7( rl) 



C F 1 



(n+l)(n + 2) 



n+l 

+4^1/j , (12) 

i=2 ' 



where Cp = 4/3. Since the moments are positive and 
increase with n, the effects of higher-order Gegenbauer 
polynomials are damped at higher scales as the DAs ap- 
proach their asymptotic form. The conformal expansion 
can thus be truncated. Quantities such as the pion's elec- 
tromagnetic form factor, for example, are given by con- 
volutions in which the kernels are slowly-varying and the 
strongly-oscillating Gegenbauer polynomials are washed 
out. The same conclusion is reached by considering, 
rather than the conformal expansion, the diagonalisation 
of the ERBL equations 0, E3, [11] which govern the evo- 



lution of the DAs much as the DGLAP equations [19H22J 
govern the evolution of PDFs. 

We can obtain values for the Gegenbauer moments 
from lattice simulations since the Gegenbauer moments 



where /+ p (0) is the vector B — s- P form factor at zero 
momentum transfer and c\ ~ 0(1). Other examples in- 
clude the ratio of the weak radiative decay amplitudes 
B — > fry and B — > K*j, where the main theoretical error 
originates from such SU(3)-breaking effects. The mea- 
sured ratio of these decay rates allows for a determination 
of the ratio of CKM matrix elements |Vf S |/|Vtd|- 

There have been three main approaches to the study 
of DAs: extraction from experimental data, calculations 
using QCD sum rules and lattice calculations. The over- 
all normalisations are given by local hadronic matrix el- 
ements, essentially the decay constants, which have al- 
ready been discussed and are partly accessible experi- 
mentally, and partly have to be calculated theoretically. 
The shapes of the leading-twist distribution amplitudes, 
in the form of the Gegenbauer moments, can be deter- 
mined from experiments by analysing data on form fac- 
tors such as F' 77 * 7r , which was studied by the CLEO ex- 
periment [25j, and the pion's electromagnetic form factor, 
F® m [26| . There is a lack of sufficiently accurate data, 
however, and it is difficult to avoid contamination from 
other hadronic uncertainties and higher twist effects. As 
a result, the existing experimental constraints are not 
very stringent. 
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Moments of DAs, then, must largely be determined 
from theory. Lattice [H, H3411| and sum rule (33l436j | 
studies have usually focussed on the second moment of 
the pion's distribution amplitude. However, the early 
lattice results were largely exploratory while sum rule 
results have an irreducible error of ~ 20% because it is 
not possible properly to isolate the hadronic states. 

The first moment of the kaon's distribution amplitude, 
for example, has in the past been determined mainly from 
QCD sum rules, and representative results include: 



of (lGeV) = 



'0.05(2) [37] 

0.10(12) [38] 

0.050(25) [39] 

0.06(3) [40] 



(16) 



These results all have the expected sign, but the un- 
certainties are around 50%. The reduction of such un- 
certainties is the chief motivation of the lattice pro- 
gramme. In an earlier publication [4ll [42[ , we obtained 
(£V(2GcV) ee 3/5<(2GeV) = 0.032(3). We note 
that in addition to the UKQCD/RBC programme for the 
calculation of DA moments on the lattice using Nf = 2+1 
domain-wall fermions, there is a UKQCD/QCDSF pro- 
gramme using Nf = 2 improved Wilson quarks [23j . 
QCDSF have also published results for moments of 
baryon DAs [28] . Lattice results for hadronic distribution 
amplitudes are considered in a recent review of hadron 
structure from lattice QCD in (43[. 

The plan for the remainder of this paper is as follows. 
In Sec. [IT] we discuss the extraction of bare moments of 
distribution amplitudes from Euclidean lattice correla- 
tion functions (we use 'bare' or 'latt' to denote quanti- 
ties before matching from the lattice to the continuum). 
In Sec. IIIII we give the details of our numerical calcula- 
tions and present the bare results. The renormalisation 
of those bare results is described in Sec. IIVI We then 
present our summary in Sec. [V] 



II. BARE MOMENTS FROM LATTICE 
CORRELATION FUNCTIONS 

In this section, we describe our general strategy for 
the lattice calculation of the unrenormalised lowest mo- 
ments of light meson distribution amplitudes. We ob- 
tain expressions for the first and second moments (f 1 ) 
and (£ 2 ) for pseudoscalar mesons and for the longitudi- 
nal moments and (£ 2 )" for vector mesons, in terms 
of Euclidean lattice correlation functions which can be 
computed by Monte Carlo integration of the QCD path 
integral. In each case, we consider a generic meson having 
valence quark content q 2 qi, where the subscripts indicate 
that the flavours of the two quarks may be different. We 
will see below that we can obtain all of these moments 
from ratios of two-point correlation functions and thus 
we expect to benefit from a significant reduction of the 
statistical fluctuations. 



A. Lattice Operators 

We now define the lattice operators used in the correla- 
tion functions from which we extract the moments of the 
distribution amplitudes. We use the following interpo- 
lating operators for the pseudoscalar and vector mesons: 



p ( x ) = 92(^)7591(2;), 

Vp{x) EE g 2 0*07/^1 0)> 

A^x) ee q^x^^q^x) . 



(17a) 
(17b) 
(17c) 



Although we have written P, V and A as local operators 
in Eq. (|17p , in the numerical simulations we use smeared 
operators at the source of our correlation functions in or- 
der to improve the overlap with the lightest meson states. 
Since the effects of smearing cancel in the ratios con- 
structed below, the discussion in this section holds for 
both smeared and local interpolating operators. We ex- 
plain the details of our smearing procedures in Scc. lIII A"l 
The operators in Eqs. © from which the moments of the 
distribution amplitudes are obtained are of course local 
operators. 

In constructing the lattice operators of Eqs. ©, we use 
the following symmetric left- and right-acting covariant 
derivatives: 

-» 1 

D fl ^p(x) = — - [U (x, x+fi)ip(x+jj) - U(x, x—p)il){x—p)\ , 



2a 



1 r- 



(18) 



ip(x)D IJy — — \Tp(x+fi)U(x+(i,x) — if>(x—jj,)U(x—jj,,x)\ , 

(19) 

where U (x, y) is the gauge link going from site x to site y 
and fl is a vector of length a in the direction /i (a denotes 
the lattice spacing). The operators of interest are then 
defined by 

<-> 

C{p^}(x) ee q 2 {x)r) {p D^qi{x) , (20a) 
C{ p ^}{x) ee q^x^yD^D^q^x) , (20b) 

°{ P n}( x ) = Q2(xh{ P l5Dp }qi (x) , (20c) 

OUW = ^)l{ P l^DpD v}qi {x) , (20d) 

where the braces in the subscripts indicate symmetri- 
sation of the enclosed Lorentz indices, {fix . . . fi n } = 

Eperms S {M S (1) • • • Ms(n)}M 



B. Operator Mixing 

In the continuum the operators in Eq. (|20l) transform 
as second- or third-rank tensors under the Lorentz group. 
On the lattice however we must consider their transfor- 
mation properties under the hypercubic group of re- 
flections and 7r/2 rotations, together with the discrete 
symmetries parity P and charge-conjugation C, where 
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the possibilities for operator mixing are increased. A 
detailed study of the transformations of these operators 
under % 4 has been performed in [44j |. 

The choice of Lorentz indices in the operators used in 
simulations is important both to keep the operator mix- 
ing simple and also to enable the extraction of matrix 
elements using as few non-zero components of momen- 
tum as possible. The latter is to avoid the associated 
discretisation effects and statistical degradation. 0{ pp } 
and C{ PAl } renormalise multiplicatively under H4 when 
p 7^ p. In the notation of [45], these operators trans- 
form under the 6-dimensional 6^ (for 0? 1 ) or 6^ _ ) 
(for 0j pp }) irreducible representations ofH.4. The choice 
p =/= p is the most convenient one for the extraction of the 
first moment of the distribution amplitudes. Charge con- 
jugation symmetry combined with "H 4 ensures that there 
is no mixing with operators containing total derivatives. 

It is also possible to obtain the first moment from the 
four operators (or 0% ■,), which each transform as 

four-dimensional reducible representations containing a 
singlet. The three traceless operators transform as the 3- 
dimensional irreducible representation (3, (without 
7 5 ) or (3, (with 7 5 ). Subtracting the trace involves 
the subtraction of a power divergence, so for the first 
moment of the distribution amplitude of the K and K* 
we avoid this by evaluating the matrix elements of 0{ pp } 
and 0{ pp } respectively with p 7^ p. 

Similarly for the second moment of the distribution 
amplitudes the most convenient choice is to use Oj 
or 0{ pM „} with all three indices different, which trans- 
form as the (1/2, l/2)(+' and (l/2,l/2)( _) 4-dimensional 
irreducible representations respectively. Charge conjuga- 



tion symmetry allows mixing of 0? pjUI/ v and 0{ PAI „} with 
operators containing total derivatives: 

°{ Pll u}{ x ) mixes with <9 {p d M (q 2 (xhu}75<li(x)) , 
0{ pliV }{x) mixes with d {p d p (q 2 (x)j u} qi(x)) . 

The moments of the distribution functions are obtained 
from non-forward matrix elements between a meson at 
non-zero four momentum and the vacuum, so the total- 
derivative operators must be included in the analysis. 



C. and (£ 2 )f from Correlation Function 

Ratios 

To obtain the first and second moments of the pseu- 
doscalar meson distribution amplitude, and (£ 2 ), we 
consider the following two-point correlation functions: 

C AvP {t,p) = Y,e tp ' x (0\A„(t,x)P\0)\0) , (21a) 

X 

C 5 M (t,p) = ^e^<0|O { 5 w} (t )a; )pt(0)|0) , (21b) 

X 
X 

For a generic pseudoscalar meson P, we define 
Zp = (P(p)\P^\0) and the bare decay constant by 
(0\A v \P{p)) ee ip„f P alc . The operators Pt(o) hiEqs. (JUJ) 
are smeared as explained below. At large Euclidean times 
t and T — t, the correlation functions defined above tend 
towards: 



CiMt p) ^ z r ^-^ llt _ T/2)Ep) ^ <22) 

C{Wt , p) (23) 



We can extract bare values for the first and second mo- 
ments of the pseudoscalar meson distribution amplitudes 
from the following ratios of correlation functions: 



R 



{p/4; 



C Au p{t,p) 



i EfiEft^i)ba«» ) (25a ) 



C Aa p(t, P ) 



PpPpVv /^2\ba 
Pa 



(25b) 



Keeping in mind the operator mixing outlined above, we 

? { P 4};4< 



obtain the first moment from R? 4 i. 4 (i,p) (the index 4 



corresponds to the time direction) with p = 1, 2 or 3 and 
a single non-zero component of momentum, \p p \ = 2tt/L. 
The second moment is extracted from R{ pll 4\.4,{t,p) with 
at least two non-zero components of momentum. We take 
p,fi = 1,2 or 3 with p 7^ p and \p p \ = \p p \ = 2n/L. We 
present more details in Sec. 1111 Bl 

Apart from isolating the moments of the DAs as much 
as possible by cancelling Zp, /b arc and most of the energy 
dependence from Eqs. (|2"3"|) and (JUJ), the ratios also sim- 
plify the effect of mixing with total derivative operators. 
These operators have matrix elements proportional to 
with which we build a ratio similar to (|25bl) . Hence 
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the contribution of the mixing term becomes trivial and 
does not have to be computed explicitly. It enters as an 
additive constant when renormalising the bare moments 
as we will discuss later. 

D. and (£ 2 )!^ from Correlation Function 

Ratios 

The treatment of the vector meson's longitudinal dis- 
tribution amplitude is analogous. We consider the fol- 
lowing two-point correlation functions: 

£7v M v-„(t,p) = XX P "<°I F m(*,*)K + (0)|0}, (26a) 

C {pp}l/ (t,p) = J2^ p - x (0\O {pp} (t,x)Vj(0)\0), (26b) 

as 

C W u}a(t,p) = Y,^^P{ P ^ } (t, X ) V ^)\Q). (26c) 



Again, the source operators V^(0) are smeared. We de- 
fine the bare longitudinal decay constant of a vector me- 
son V , with polarisation index A and polarisation vector 
ei X \ by (0|V„|V(p,A)) = f v ^ c m v e {X) . Then, at large 
Euclidean times t and T — t, the correlation functions 
defined above may be written: 



r (tn) , -(fr c ™v) 2 z- EvT/2 ™M(t-T/2)E v ) ( p, Pu \ 

^VpvAhP) _ ~9nu H — ' \ 21& ) 

E v \ m v J 

-»(/barc my) 2 e - gv T/2 (g l ) || bare^^ T/2 )E V ) 1 / 2p P Pp,Pu \ 

C{ pp } v (t,p) -> — ^-^Pp - g^Pp + —-3 — j , (27b) 

(/b-mv) 2 e ^^ T / 2 (g 2 )H barc sinh((t - T/2)£ y ) 1 / ^ 3 WPi/Po 

(27c) 

where we have used the completeness relation for the polarisation vectors of massive vector particles, ^\ e^e^^ = 
9^1!^ ~T" Pp.Pl/ /mfr. We extract bare values for the first and second moments from the following ratios: 



s -> -i<^>ll — tanh((i - T/2)E V ) \ (-g pu p, - g, uPp + 2 -MfA 



, (28a) 



3 

nV _\ = C {p p v}a (t,p) 

{P ^ }A ~ lY.iCv^Vi = o,\p\ = ¥ ) 



-> (£ 2 ) 11 barc tanh((t - T/2)E V ) \ (-g paPpPu - g^p p p„ - g vaPpPp + ^ ) , ! 2sl , , 



where the index i runs over spatial dimensions only. We 
obtain the first moment from R\ p Q v {t,p) at p = by 
taking p = ^ = l,2or3. The second moment is obtained 
from R{ ppL ,} a (t,P) by taking, for example, v = 4,p = 
1 , /i = a = 2 and a single non-zero component of p in the 
1-direction. 



III. NUMERICAL SIMULATIONS AND 
RESULTS 

A. Simulation Details 

Our numerical calculations are based upon gauge field 
configurations drawn from the joint datasets used for 
the broader UKQCD/RBC domain- wall fermion phe- 
nomenology programme. Configurations were gener- 
ated with Nf = 2 + 1 flavours of dynamical domain- 
wall fermions and with the Iwasaki gauge action, us- 
ing the Rational Hybrid Monte Carlo (RHMC) [t| 



7 



algorithm on QCDOC computers [47H49] running the 
Columbia Physics System (CPS) software [50j and the 
BAGEL [5l], [52j assembler generator. 

Our set of gauge configurations includes data with two 
different volumes but at a single lattice spacing, thus giv- 
ing us some indication of the size of finite volume effects 
but no ability to perform a continuum extrapolation. We 
therefore have an unavoidable systematic error which is, 
however, formally of 0(a 2 Aq CD ) due to the automatic 
O(o)-improvement of the DWF action and operators. In 
the future, this limitation will be overcome by perform- 
ing the analysis with a dataset with a finer lattice spac- 
ing with the same action (such a dataset is now available 
and is currently being calibrated). In the meantime, fol- 
lowing the UKQCD/RBC procedure for these configura- 
tions we ascribe a 4% uncertainty as the discretisa- 
tion error on the moments. For both volumes, we have 
a single dynamical strange quark mass, close to its phys- 
ical value. We use several independent ensembles with 
different light-quark masses (m u — rrid), all heavier than 
those found in nature. The hadronic spectrum and other 
properties of these configurations have been studied in 
detail and the results have been presented in [53J (for the 
lattice volume (L/a) 3 x T/a = 16 3 x 32) and O (for the 
lattice volume 24 3 x 64) . In both cases the length of the 
fifth dimension is L s = 16. 

The choice of bare parameters in our simulations is 
P = 2.13 for the bare gauge coupling, am s = 0.04 for the 
strange quark mass and am q = 0.03, 0.02, 0.01 and, in 
the 24 3 case only, 0.005 for the bare light-quark masses. 
A posteriori, the strange quark mass is found to be about 
15% larger than its physical value. The lattice spacing 
is found to be a" 1 = 1.729(28) GeV [H, giving physi- 
cal volumes of (1.83 fm) 3 and (2.74 fm) 3 . The lattice 
spacing and physical quark masses were obtained using 
the masses of the 7r and K pseudoscalar mesons and the 
triply-strange Q baryon. The quark masses obtained in 
the 24 3 study are shown in Table [J] Owing to the rem- 
nant chiral symmetry breaking, the quark mass has to 
be corrected additively by the residual mass in the chiral 
limit, am xes = 0.00315(2) [l3|]. The physical pion masses 
are as follows: 



670 MeV am 
555 MeV 
415 MeV 



0.03 
am q = 0.02 
am q = 0.01 
330 MeV am q = 0.005 



(29) 



differs from that used in [53j in that the Markov chains 
have been extended for the heaviest light quark mass to 
give additional statistics, using an improved algorithm 
that decorrelated topology rather more quickly. 

In order to improve the statistical sampling of the cor- 
relation functions, on each configuration we have aver- 
aged the results obtained from either 2, 3 or 4 sources 
spaced out along a lattice diagonal. In the 16 3 case, for 
example, the sources used are at the origin, at (4, 4, 4, 8), 
(8, 8, 8, 16) and (12, 12, 12, 24). Statistical errors for ob- 
servables are estimated using single-elimination jack- 
knife, with measurements made on the same configura- 
tion but at different source positions put in the same 
jackknife bin because of the correlations expected be- 
tween them. In order to lessen the effect of autocorrela- 
tions, we follow the same blocking procedures as in (53j 
and [13. In the 16 3 case, the span of the measurements 
in each block covers 50 molecular dynamics time units. 
In the 24 3 case, for the m q a = 0.005 and am q = 0.01 en- 
sembles, each jackknife bin contains measurements from 
every 80 molecular dynamics time units, while for the 
am q = 0.02 and am q — 0.03 ensembles each bin contains 
measurements from every 40 molecular dynamics time 
units in order to have a reasonable number of bins for 
the analysis. 

We use source smearing to improve the overlap with 
the mesonic states, either gauge-fixed hydrogen S- 
wavefunction smearing [55[ with radius r — 3.5 in lat- 
tice units or gauge invariant Gaussian smearing |56j with 
radius r — 4. 



B. Results 



In order to extract (S, 1 )k from the ratio R{ pfJi y v (t,p) 



defined in (|25a[) . we need the two correlation functions to 
be measured at \p\ 7^ 0. Since we expect hadronic observ- 
ables with larger lattice momenta to have larger lattice 
artefacts and statistical errors, we restrict the choice of 
indices to p = v = 4 and fi = 1, 2 or 3 with \p\ = 2ir/L 
(i.e., p M = ±2-7r / L with the remaining two components 
of p equal to 0). (S, 1 )^ 10 can then be obtained from the 
ratio at large times: 



2 k 



flf«};4(*>Pfc=±27r/£) = ±i-H£ 



1 \ bare 



(30) 



Measurements were performed using the UKhadron 
software package that makes use of both the BAGEL 
DWF inverter [H [H[ and elements of the SciDAC 
software library stack including the Chroma LQCD li- 
brary [54| and QDPH — K The details are summarised in 
Tables [TT| and IIIII We restrict our analysis to the uni- 
tary data for which the valence and sea quark masses are 
the same (partially-quenched data was used extensively 
in the studies of the chiral behaviour of the spectrum and 
decay constants in [13]). On the 16 3 lattice, our dataset 



with \p\ = 2-k/L and k = 1, 2, 3. The plots in Fig.[T]show 
our results for (^ 1 )^ ire as a function of t obtained from 
the ratio Rf ik y ±{t,Pk= ± 2n/L) for the four values of 
the light-quark mass, combining results at t with those 
at T — t — 1. The results have been averaged over the 
three values of k and, in total, the 6 equivalent lattice 
momenta with \p\ = 2n/L. 

To obtain (C 2 )^ f™ the ratio R? ptlv} . a {t,p) defined 
in (I25b[) we need two non-zero components of momentum, 
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TABLE I. Lattice scale and unrenormalised quark masses in lattice units, from the 24 4 lattices [Ti§| • Note fhx = mx + m t , 
Only the statistical errors are given here. 



a' 1 [GeV] 


a [fm] 


am ud 


am ud 


am s 


am s 


arrivd ■ am s 


1.729(28) 


0.1141(18) 


-0.001847(58) 


0.001300(58) 


0.0343(16) 


0.0375(16) 


1:28.8(4) 



am q 
am q 



= 0.01 h~h 
= 0.02 «h 
= 0.03 * 



1(1 



14 



l(i 



0.04 
0.035 

0.03 
0.025 

0.02 
0.015 

0.01 



am, = 0.005 
am q — 0.01 
am q — 0.02 
am a — 0.03 



0.005 f „ * * 1 * 1 * $ 




10 



15 



211 



25 



30 



FIG. 1. Results for (C 1 )^ 10 as a function of the time, on the 16 3 (left) and 24 3 (right) lattices. The shaded band shows the fit 
range, fitted value and its error. 



TABLE II. Parameters for our 16 3 dataset, which corresponds 
largely to that of [53| . The range and measurement separation 
A are specified in molecular dynamics time units. Nm eas is the 
number of measurements for each source position f src . The to- 
tal number of measurements is therefore iV mcas x N alc , where 
TVsrc is the number of different values for t src . In the right- 
most column, XY-XY denotes contraction of two quark prop- 
agators with X-type smearing at source and Y-type smearing 
at sink: G = Gaussian wavefunction, L = point. 



mi 


Range 


A 


me as 


tsrc locations 


Smearing 


0.01 


500-3990 


10 


350 


0,8, 16,24 


GL-GL 


0.02 


500-3990 


10 


350 


0,8, 16,24 


GL-GL 


0.03 


4030-7600 


10 


358 


0,16 


GL-GL 



TABLE III. Parameters for our 24 3 dataset, which corre- 
sponds to the unitary part of the dataset of 13]. Columns 
as in Table [TT] with addition of H = gauge-fixed hydrogen S- 
wave smearing. 



mi 


Range 


A 


■^v me as 


t S rc locations 


Smearing 


0.005 


900-4480 


20 


180 


0,32,16 


HL-HL 


0.01 


800-3940 


10 


315 


0,32 


GL-GL 


0.02 


1800-3580 


20 


90 


0,32 


HL-HL 


0.03 


1260-3040 


20 


90 


0,32 


HL-HL 



so we use 

R{4jk};4(t,Pj = ±^/L,p k = ±2n/L) = 

-(±^)(±^)<a bare , (31) 



with \p\ — ^/22tt/L, k,j — 1, 2, 3 and k ^ j. We average 
over all 4 momentum combinations appropriate to each 
of the 3 possible Lorentz index choices. 

We may extract (£ 1 )l : » arc from the ratio R^ p ^ u (t,p) 
defined in (|28aj) by considering only zero-momentum 
correlation functions. In the denominator, we average 
CviVi {t, P = 0) over all 3 spatial directions. In the numer- 
ator, we average over C{ 41 } 1 (t,p = 0), C{ 42 }2(^P = 0) 
and C{4 3 } 3 (t,p = 0). Results are shown in Fig. [5] 

(£ 2 )|-* a p C i s extracted from the ratio defined in (|28bl) by 
averaging C^v«(t,Pi = 0, \p\ = ^ ) over all 4 appropri- 
ate momenta for all 3 spatial directions in the denomina- 
tor. In the numerator we average over all possible com- 
binations of C^ijyi(t,pj = \p\ = with i j. In 
principle we should include disconnected contributions in 
the (f> correlation functions. We argue that these contri- 
butions are Zweig-suppressed however and can therefore 
be neglected. 

If picking the fit range was not straightforward, we 
considered the correlation functions in the numerator 
and denominator separately. We identified and excluded 
from our fits the region where the excited states still con- 
tributed. We then chose the fit range aiming for a good 
X 2 /d.o.f. and a stable fit with respect to small variations 
of the lower bound of the range. Owing to the increasing 
noise when t gets larger, the fits are insensitive to the 
upper bound of the fit range. 

The 16 3 and 24 3 bare results are given in Tables HVl and 
|V] respectively, complete with linear chiral extrapolations 
which, as we shall discuss in the next section, can be 
justified using chiral perturbation theory (at least in the 
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0.04 
0.035 

0.03 
0.025 

0.02 
0.015 

0.01 
0.005 




I 

i I 



15 



211 



25 



30 



FIG. 2. Results for (£ )| f » arG as a function of the time, on the 16 (left) and 24 (right) lattices. Symbols as in Fig. [T] 



pseudoscalar case). 



C. Quark Mass Extrapolations 

In leading-order chiral perturbation theory [57|, 
is proportional to m s — m u u without chiral logarithms: 



(a 1 : 



K 



8B 
P 



j,/d)h 



(32) 



where / and Bq denote the usual chiral perturbation the- 
ory parameters and 6i : 2 is a Wilson coefficient introduced 
in [53. Our data shows clearly the effects of SU(3) sym- 
metry breaking and is compatible with this expectation. 
We therefore perform a linear extrapolation in a(m s — m q ) 
to the physical point a(m s — m u( j), as shown in Fig. [3] 
The second error quoted in the results in the chiral-limit 
for the first moments in Tables IIVI and [V] is due to the 
uncertainty in this physical point (determined using the 
quark masses in Table [TJ) . In this way we deal simul- 
taneously with the usual light-quark mass extrapolation 
and with the strange quark mass extrapolation which is 
necessitated by our strange quark mass being approxi- 
mately 15% too heavy. We have not constrained our fit 
to vanish in the SU(3) limit. 

A similar linear behaviour is seen for (£ 1 )1 f » aro (see 
Fig. [3]), so we follow the same extrapolation procedure. 
We note a hint of a finite volume effect in the K* case but 
not in the K case, which is contrary to what we would 
expect. Where we have K* results for both volumes at 
the same light-quark mass, however, they agree within 
the statistical uncertainties. 

For the second moments, we also have some guidance 
from chiral perturbation theory [58) : there is no non- 
analytic dependence at 1 loop and we should fit linearly 
in m^.. The dependence on the quark masses is very mild 
in any case and in fact our results for the p, K* and <p 
agree within the statistical errors. Therefore we perform 



a linear extrapolation in the light quark masses and ne- 
glect the effect of the too-heavy strange quark mass (see 
Fig. [4]). We see no indication for finite size effects in 
the second moments when we compare the data points 
on the two different lattice volumes. They agree within 
their statistical errors. 



IV. RENORMALISATION OF THE LATTICE 
COMPOSITE OPERATORS 

We now discuss the conversion of our bare lattice re- 
sults to results in the MS scheme. To reduce system- 
atic uncertainties we have determined the renormalisa- 
tion factors nonperturbatively in the RI'/MOM scheme, 
continuing the work in [59|, and convert to MS using 3- 
loop continuum perturbation theory [6fJ, l61| . We begin, 
however, with a perturbative calculation of the renor- 
malisation factors. T he p erturbative results have been 
used previously in [62, 63] and will provide a comparison 
to the nonperturbative results. The contribution to the 
second moment from mixing with a total-derivative op- 
erator is calculated perturbatively only. We will see that 
this contribution is small and is not accessible within the 
current nonperturbative scheme. 



A. Perturbative Renormalisation 

The perturbative matching from the lattice to MS 
schemes is performed by comparing one-loop calculations 
of quark two-point one particle irreducible (1PI) func- 
tions with an insertion of the relevant bare lattice opera- 
tor. This requires the evaluation of the diagrams shown 
in Fig. [SI together with wavefunction renormalisation fac- 
tors, Fig. [5] For the first-moment operator, we define 



(33) 



10 




0.035 



0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 

am s — am q 




0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 

ani g — am,. 



FIG. 3. Chiral extrapolations for (£ )k ic an d {£ )it* arC - The extrapolation to the physical point is shown by the vertical solid 
line, with uncertainty, dominated by the uncertainty in the physical strange mass, indicated by the dotted lines. 



TABLE IV. Summary of results for the bare values of the distribution amplitude moments on the 16 lattices. The chiral 
extrapolations are discussed in Sec. MI CI and the errors are statistical and (in the first moment case) due to the uncertainty 
in the physical point for the chiral extrapolation. 



am ud 


0.03 


0.02 


0.01 


0.005 


X-limit 




0.110(2) 


0.109(2) 


0.113(4) 




0.112(5) 


/ /-l \ bare 
\? IK 


0.00543(27) 


0.01174(71) 


0.0194(15) 




0.0228(14)(11) 


1 (-2\ bare 
\? IK 


0.109(2) 


0.107(2) 


0.113(3) 




0.112(4) 


(c 2 >I barc 


0.113(4) 


0.100(5) 


0.116(6) 




0.109(10) 




0.00610(24) 


0.01275(51) 


0.0207(10) 




0.02443(96) (107) 




0.111(4) 


0.101(4) 


0.113(4) 




0.110(6) 


<e 2 >r arc 


0.109(3) 


0.100(3) 


0.109(3) 




0.107(5) 



For the second moment calculation we must take account 



of mixing with a total derivative operator (c.f. Sec. IIIBj) 
Adopting the notation 



Odd = ^{f/rsDuD^tp, O dd = d^d^^^, (34) 



with all Lorentz indices distinct and symmetrised, we 
need to determine 



OdSM = Z D D,DD^a)O l ^{a) + Z D D.dd(HOod (a)- 
The renormalisation factors are given by 



Z ^M)= {1 _ wl)Zw 

ZDD.DD^a) = {1 _ w2q)Zw 

ZddmA^o) = 1 2 [\ Mpa) + V c ^ s - Vas 

(l — W )Z w 47T \6 



aC-p 
Air 
aC-p 

47T 



1G 

T 

25 



ln(/xa) + Ef s - Si + V MS - V 



-— ln(/m) + £f b - Ei + Vgg - V DD 



(36) 
(37) 
(38) 



In the equations above (1 — Wq)Z w is a characteristic form at one-loop as 
normalisation factor for the physical quark fields in the 
domain- wall formalism. Z w represents an additive renor- 
malisation of the large Dirac mass or domain- wall height aCV 2wq 
AI = 1 - w , which can be rewritten in multiplicative Z w = 1 + z w , z w = - - ^ 2 T, w . (39) 
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0.1!") 



0.15 



16 J x 32 i — A- 
24 3 x 64 i — 





0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 
arriq + am res 



0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 
am q + uriires 





0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 
am q + am r es 
0.15 



0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 
am q + am Ies 



0.14 - 
0.13 
0.12 
0.11 
0.1 
0.09 
0.08 
0.07 - 
0.06 



16 3 x 32 i — A- 
24 :i x 64 i — t- 




0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 
am q + am res 



FIG. 4. Chiral extrapolations for (£ 2 )° arc , {£, 2 }k, 
shown by the solid vertical line in each case. 



c -2\ || bare / >-2\ || bare 



and (£ z ) 



2\ || bare 



The physical value for am q + am ICS is 



The one-loop correction z w becomes very large for cer- 
tain choices of M [64], HI], including that used in our 
numerical simulations, so that some form of mean-field 
improvement is necessary, as discussed below. 

Terms with superscripts MS in Eqs. ([36]), <[37j) and ([38]) 
arise from the continuum calculations, while unsuper- 



scripted terms come from the computations in the lat- 
tice scheme. To shorten some expressions below we will 
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TABLE V. Summary of results for the bare values of the distribution amplitude moments on the 24^ lattices 



am ud 


0.03 


0.02 


0.01 


0.005 


X-limit 


I c'l \ bare 
\S lit 


n i n^i'Q'i 

u. luo^y j 










/ £l \ bare 
\? IK 


0.00566(33) 


0.01254(72) 


0.01946(65) 


0.0231(15) 


0.02377(71)(110) 


1 (-2\ bare 
\4 IK 


0.103(8) 


0.106(4) 


0.112(2) 


0.113(6) 


0.117(5) 


(c 2 )H barc 


0.110(9) 


0.093(10) 


0.112(3) 


0.120(13) 


0.118(7) 




0.00619(35) 


0.0139(10) 


0.0225(13) 


0.0311(30) 


0.0281(13)(14) 




0.109(12) 


0.095(8) 


0.108(3) 


0.117(5) 


0.118(7) 


(c2) ybare 


0.108(7) 


0.097(7) 


0.105(2) 


0.107(3) 


0.107(4) 



TABLE VI. Constants needed for the perturbative renormalisation of the first and second moment operators using domain-wall 
fermions and the Iwasaki gauge action (ci = —0.331). M is the domain-wall height, c = E^ s — Ei + V MS — V, cdd = 
yjMS _ _|_ Vdd — Vdd and egg = Vgg S — Vgg. Ei, V , Vdd and Vgg are dependent on the gauge and the infrared regulator: 
Feynman gauge and a gluon mass are used here. V was calculated in [4l[, while Vdd and Vgg have been calculated as part of 
this work. 



M 


Ei 


V 


c 


Vdd 


CDD 


Vgg 


egg 


0.1 


4.6519 


-4.6297 


-0.9110 


-10.816 


4.9838 


0.5415 


0.0279 


0.2 


4.5193 


-4.5614 


-0.8468 


-10.698 


4.9982 


0.4285 


0.1409 


0.3 


4.4093 


-4.5101 


-0.7881 


-10.608 


5.0179 


0.3433 


0.2262 


0.4 


4.3158 


-4.4678 


-0.7369 


-10.533 


5.0362 


0.2729 


0.2966 


0.5 


4.2354 


-4.4311 


-0.6932 


-10.467 


5.0509 


0.2119 


0.3575 


0.6 


4.1665 


-4.3980 


-0.6574 


-10.407 


5.0603 


0.1573 


0.4122 


0.7 


4.1079 


-4.3673 


-0.6295 


-10.352 


5.0639 


0.1070 


0.4625 


0.8 


4.0593 


-4.3381 


-0.6101 


-10.300 


5.0604 


0.0597 


0.5098 


0.9 


4.0204 


-4.3097 


-0.5996 


-10.250 


5.0489 


0.0142 


0.5552 


1. 


3.9915 


-4.2816 


-0.5988 


-10.200 


5.0283 


-0.0303 


0.5998 


1.1 


3.9731 


-4.2529 


-0.6090 


-10.151 


4.9970 


-0.0749 


0.6443 


1.2 


3.9664 


-4.2232 


-0.6321 


-10.100 


4.9528 


-0.1205 


0.6899 


1.3 


3.9727 


-4.1916 


-0.6700 


-10.047 


4.8933 


-0.1682 


0.7376 


1.4 


3.9943 


-4.1571 


-0.7261 


-9.9895 


4.8147 


-0.2195 


0.7889 


1.5 


4.0343 


-4.1182 


-0.8050 


-9.9267 


4.7119 


-0.2764 


0.8458 


1.6 


4.0974 


-4.0728 


-0.9135 


-9.8551 


4.5771 


-0.3418 


0.9112 


1.7 


4.1905 


-4.0176 


-1.0618 


-9.7700 


4.3989 


-0.4205 


0.9899 


1.8 


4.3249 


-3.9462 


-1.2676 


-9.6627 


4.1572 


-0.5211 


1.0905 


1.9 


4.5209 


-3.8447 


-1.5651 


-9.5140 


3.8125 


-0.6631 


1.2325 



define 

c=zf E -Z 1 + V m -V, (40) 

con = zY S -Zi + V™E-V D D, (41) 

c dd = V^-V ad . (42) 

The terms s and Si come from quark wavefunction 
renormalisation, while V MS , Vfffi, V^ s and V, Vdd, 
Vqq come from the one-loop corrections to the ampu- 
tated two-point function. They are given by "vertex" 
and "sail" diagrams, plus an operator tadpole diagram 
in the lattice case. V^fp and Vdd can be isolated by 
computing the one-loop correction with equal incoming 
and outgoing quark momenta. Likewise Vqq S and Vqq are 
found by setting the incoming and outgoing quark mo- 
menta equal and opposite (the lattice tadpole diagram 
does not contribute in this case). Using naive dimen- 
sional regularisation (NDR) in Feynman gauge with a 




FIG. 5. One-loop vertex diagrams evaluated in the perturba- 
tive renormalisation of the 1st and 2nd moment operators. 




FIG. 6. One-loop diagrams for the quarks' wavefunction 
renormalisation. 
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TABLE VII. Values for z w , 
in [6^|, and dj extracted from [7l| 



z;„" extracted from the results 



M 




MF 




0.1 


—243.86 


—86.579 


—0 09303 


0.2 


— 113.29 


—39.501 


—0 01 7Q8 


0.3 


—69.404 


—93 830 


—0 014Q7 


0.4 


—47.077 


— 15.949 


—0.01274 


0.5 


—33.278 


— 11.142 


—0 01 0Q0 

VJ . VJ lUt/U 


0.6 


—23.648 


—7.8365 


—0 00Q31 5 

U . VJVJUO -LO 


0.7 


— 16.300 


— 5 3538 


—0 0078Q.fi 


0.8 


— 10.263 


—3.3459 


—0 00658Q 


0.9 


—4.9617 


— 1.6078 


—0 005379 

vj • vjvjtj'J I iJ 


1.0 


0.0 


0.0 


—0 004961 

VJ . VJVJ'rtZiVJ _L 


1.1 


4.9442 


1.5902 


-0.003227 


1.2 


10.192 


3.2748 


-0.002290 


1.3 


16.136 


5.1900 


-0.001485 


1.4 


23.346 


7.5350 


-0.0008650 


1.5 


32.784 


10.648 


-0.0005360 


1.6 


46.322 


15.194 


-0.0006566 


1.7 


68.294 


22.720 


-0.001570 


1.0 


111.03 


6i .yui 


— U.UU4U14 


1.9 


241.55 


84.270 


-0.01020 


gluon 


mass infrared (IR) rej 


;ulator, 






yMS 1 

1 2' 


yMS _ 


<«> 




T/ MS 121 

Vdd- 


t/MS 





The lattice contributions are evaluated for domain-wall 
fcrmions with the Iwasaki gluon action (ci = —0.331), 
also choosing Feynman gauge and using a gluon mass 
IR regulator. Ei has been evaluated in 65], while we 
calculated the vertex term V for the first moment op- 
erator in [4l| . Here we have calculated the vertex terms 
Vd d and Vgg for the second moment operator. Perturba- 
tive calculations with domain- wall fermions are explained 
in [65j and the form of the Iwasaki gluon propagator 
can be found in [66]. Values for Si, V, Vdd and Vgg are 
given as functions of M in Table fVTl along with c, cdd 
and egg. Chiral symmetry of the domain- wall action im- 
plies that these results also apply for the operators which 

I 



are like those used here, but without the 75. We note that 
the perturbative renormalisation factor for the first mo- 
ment operator using alternative fermion and gauge for- 
mulations can be found in [67i ] (domain- wall fermions and 
plaquette action), l68l (overlap fermions and Liischer- 
Weisz action) and [69] (clover fermions and plaquette 
action). Second moment calculations with clover and 
Wilson fermions have been performed in [6!| and [7(3] 
respectively (in both cases using the plaquette action). 

Our numerical simulations use M = 1.8. For this value 
of M, with the Iwasaki gluon action, the one-loop coef- 
ficient in the physical quark normalisation is z w ~ 112 
(extracted from T, w in Table III of (65[), making it clear 
that mean-field improvement is necessary. We follow the 
prescription given in [651 ]. The first step is to define a 
mean-field value for the domain-wall height, 



M 



MF 



M - 4(1 - P 1/4 ) = 1.3029 



(45) 



where P = 0.58813(4) is the average plaquette value in 
the chiral limit in our simulations. The physical quark 
normalisation factor becomes [l — (w^ F ) 2 ] Z^ F , with 



/•MF 



.MF 



47T 

2w™ F 



MF 



1 - {w™ F ) 2 



(E w + 32n 2 T MF ) = 5.2509, 



(46) 



0.0525664 65] is a mean-field tadpole fac- 



where Tmf 

tor and T, w is evaluated at M MF . Values for z^ F as a 
function of M are quoted in Table IVII1 extracted from 
the results in [gj|. Likewise, £1 = 3.9731, V = -4.1907, 
Vdd = —10.045 and Vgg — —0.1696 are evaluated at 
A/ MF . 

For the operator Odd with two covariant derivatives, 
mean- field improvement introduces a factor ii pt / u where 
u is the mean link (here taken to be u = P 1 / 4 ) and 



1 



MF 



is its perturbative expansion. For Ogg with two ordinary 
derivatives, in contrast, the extra factor is u/u p t- The 
mean-field-improved matching factors are thus 



Z, 



MF 



z 



MF 

DD.DD 



7MF 



DD,dd 



1 



K MF ) 2 
1 1 



1 
1 

U 1 - (w™ F ) 2 ZMF 
1 



47T 

+ —— 

4ir 



10 



ln(/ia) 



MF 



25 

._ln(/xa) + c^-87r 2 r M F 



1-( W ^ F ) 2 Z™ F 4 ^ ^ 3 



1 aC F ( 5 mf 

- ln(/Xa) + Cgg 



(47) 
(48) 
(49) 



r 



with 

JV1F . 
c dd 



„MF 



,MF 



c— = -0.6713, c DD 



87t 2 Tmf = 0.7408 and two choices for the coupling. The first is a mean-field 



0.7391. To evaluate these expressions, we make 
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improved coupling defined using the measured plaquette 
value P, according to [Til] 



1 



5mf(m) 5o 



P , 
-a +d„ 



22 



ln(//a) 
4 



48tt 2 



ln(/z<z) 



(50) 



where Nf is the number of dynamical quark flavours. For 
the Iwasaki gauge action with c\ = —0.331, the values 
d g = 0.1053 and c p = 0.1401 are given in (65[, while 
values for df as a function of M were calculated in [7l[ 
and are quoted in Table I VIII In our simulations, j3 = 
6/g 2 = 2.13 with N f = 3 and a^ 1 = 1.729 GeV. The 
second choice is the continuum MS coupling, calculated 
as outlined in Appendix A of (59[. At fia — 1, we find 

a MF = 0.1769 and a m = 0.3138. We use these two 
values to evaluate the renormalisation factors above. We 
also evaluate the mean-field improved expression for the 
axial vector current renormalisation [65] , interpolating to 
our mean-field Af MF . The values are shown in Table IViTTI 
The ratios of the renormalisation factors, from which the 
factor 1/(1 — (wq) 2 )Z^ cancels, are also shown in the 
table. 

We take the mean value of the results with the two 
different choices for the coupling as the best answer for 
the renormalisation factors. The difference between the 
two choices will form the error. The relevant factors for 
the perturbative renormalisation of the ratios in Eqs. (|25[) 
and ((28]) are given in Table HXl Chiral symmetry here en- 
sures that we do not have to distinguish between vector 
and axial- vector operators. We note that the contribu- 
tion from the mixing term Zdd.qq is smaller than the 
error on Zdd dd itself. 



B. Nonperturbative Renormalisation 

In order to renormalise the correlation functions non- 
perturbatively we make use of the Rome-Southampton 
RI'/MOM scheme [72j which we now briefly review and 
discuss refinements to [59j . The starting point and def- 
inition of the RI'/MOM scheme is a simple renormali- 
sation condition that can be imposed independently of 
the regularisation used, thus on the lattice as well as in 
the continuum. This facilitates scheme changes which is 
important for the matching to MS. The renormalisation 
condition has the form 



RI'/MOM scheme these momenta are exceptional, that 
is equal incoming and outgoing quark momenta, p and 
p' . For some renormalisation factors it is advantageous 
to use a non-exceptional symmetric choice of momenta 
p 2 = p' 2 = q 2 , where q = p — p' , leading to the dis- 
tinct RI/SMOM scheme. This suppresses unwanted in- 
frared effects in the vertex amplitude, pion poles for ex- 
ample, and suggests a better-behaved accompanying con- 
tinuum perturbation theory 73]. Exceptional momenta 
with q = also cause matrix elements of operators with 
total derivatives to vanish, making Zr)D,dd inaccessible 
in our nonperturbative analysis. 

The vertex amplitude is constructed from the unam- 
putated Green's function 



Go(p) = (iP(p)O(0)i>(p)) , 
O(0) =J2^)M^x')i/j(x') 



(52) 



The external quark lines need gauge fixing, for which we 
use Landau gauge. The current J has the appropriate 
Dirac structure and may be non-local if the operator con- 
tains derivatives. For example, a single right derivative 

— y 

D v in the vector case would correspond to 

J o Pfl (x,x') = j p -(u(x,x')6 X ' >x+ p, - U(x,x')6 x ' tX -fij 

(53) 

matching the definition in Eq. (Tl"8"l) . 

The vertex amplitude itself is found after amputating 
the Green's function and tracing with a suitable projector 



A (p) = Tr[IL (p)P } (54) 



with 



U (p) = (Sip))- 1 (G (p)) (Sip))- 1 . (55) 

We have used the quark propagator S(p) and the an- 
gle brackets indicate the gauge average. The projector 
Pq depends on the particular operator and includes an 
overall normalisation factor to account for the colour and 
Dirac trace. In a simple example Pq would isolate the 
tree-level contribution to the vertex amplitude; we will 
detail our choices below. We have now defined the renor- 
malisation procedure and will turn to details of the im- 
plementation before discussing the results. 



Momentum sources 



A 



o(jp) = Z {n)Z-\v)k h ™ 



(P) 



Ag c °(p), (51) 



where Aq (A^f IC ) is the renormalised (bare) vertex am- 
plitude. Together with the quark field renormalisation 



Z q , defined by ip = Z q -0 baro , this defines the renor- 
malisation constant Zq for the operator O. The renor- 
malisation scale fj, is set by the momentum of the exter- 
nal states entering the vertex amplitude. In the original 



One refinement to our previous work [59| is the use of 
momentum sources [74j . In contrast to the point sources 
used before, this effectively amounts to a volume aver- 
age over the lattice resulting in much smaller statistical 
errors [63]. Starting from (|5"2")l the Green's function in 
momentum space is 

G o(p) = E (75^b),75 Jo{x,x') S{p) x ,) , (56) 
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TABLE VIII. Perturbative renormalisation factors and their ratios for two choices of the strong coupling, evaluated at pa = 1. 





7 MF 


7 MF 
^DD,DD 


7 MF 


7 MF 
Z'A 


Z MF 
°{pp} 
Z A 


7 MF 
Z A 


7 MF 
7 M1.' 


IMF 


0.9896 


1.1604 


0.0122 


0.8009 


1.2356 


1.4488 


0.0152 


a™ 


0.9162 


1.0966 


0.0202 


0.6934 


1.3214 


1.5815 


0.0291 



TABLE IX. Perturbative renormalisation factors to match the 
lattice results to MS at ap = 1. 



z A 



Z DD,DD 
Za 



Za 



1.28(4) 1.52(7) 0.022(7) 



where rather than use the quark propagator S(x\y) ob- 
tained by inverting the Dirac Matrix M on a point source 



Y J M{x\x)S(x\y)=5 x 



',y ■ 



(57) 



we use S(p) x — S(x\y)e lpy which can be found by 
inverting with a momentum source [zij 



Mix 1 »5(p). 



(58) 



and is defined on all lattice sites corresponding to the 
off-shell quarks used in the Green's function. The gain in 
statistical accuracy is paid for with a separate inversion 
for every momentum used in the simulation. However, 
this is more than compensated by a much reduced num- 
ber of necessary configurations. Limiting ourselves to a 
few carefully chosen momenta, statistical fluctuations are 
reduced with lower overall computational cost. 

The momenta we use are first of all constrained to be 
within a range A QC d < P 2 < 1/a for the RI'/MOM 
scheme 72 1. We use our previous results 59] to identify 
suitable values and focus on momenta which are expected 
to have reduced hypercubic lattice artefacts by trying to 
limit YlPu, f° r fixed p 2 [63[ (see also (zEIlll). The values 
used are: 

16 3 x 32: (1,1,2,3), (1,1,2,4), (1,2,2,1), 
(1,2,2,3), (1,2,2,4) 

and 

24 3 x 64 : (2,2,2,7), (2,2,2,8), (2,2,3,7), 
(2,2,3,8), (2,3,3,7), 

where we have given for momenta p^ = ln-ri^fL (with 
L — > T for time components). 



2. Projectors 

We extend the set of operators considered previously 
in [59] . We now require operators with up to two deriva- 



tives, Or^ j (n < 3), making the the necessary pro- 
jectors Pq slightly more involved than for bilinears. Since 
we resort to readily available calculations [6(J 0, IZZl 
for the final conversion to MS as well as to account 
for running, we have to tailor the projectors to match 
the RI'/MOM scheme and vertex functions used in the 
continuum calculations. Decomposing the amputated 
Green's function into terms allowed by Lorentz symme- 
try and rememberin g th at we are taking all indices to be 
distinct, we find IsaM Inl 

G a (p) = Ei(p)7 {Atl p M2 . . .p^y + T, 2 (p)p^ . . -P^tf- 

(59) 

For simplicity we limit the discussion to the vector 
case here; axial-vector operators are analogous. The 
RI'/MOM scheme uses the contribution from only 
in (|59")) . The required projector Pq will depend on the 
momentum entering the Green's function and its (fixed) 
directions fa (i = 1 . . . n). In general, multiplying Go 
with 7^ picks up combinations of both terms Ei and £2. 
On the other hand, projecting with j p where p <£ {/J-i} 
is only sensitive to £2 (note that we have n < 3). Thus 
multiplying with the difference of the two Dirac matrices 
with appropriate normalisation and momentum factors 
ensures that the vertex amplitude in (|54l) contains £r(p) 
only. There are simpler special cases in which one or more 
components of the momentum p are zero, causing the sec- 
ond term in (I59[) to vanish. However, since we tried to 
choose our momentum directions close to the diagonal of 
the lattice, we do not have momentum components that 
are zero. 

For fixed indices fii (i = 1 . . . n) of the Green's function 
we can construct n different projectors Po,i by starting 
from any of the 7^ : 

Put 

7 W - 7 P — 

Po,i= n Vp , with i = 1 . . . n . (60) 

^ n ph 

The normalisation AT is chosen such that for the tree-level 
vertex amplitude we find A^ cc (p) — 1 . The index p is dif- 
ferent from any of the p,i and such that its momentum 
component p p is as small as possible to reduce discreti- 
sation errors. We use = sinp^ to better account for 
lattice momenta. The case of axial-vector operators O 5 
is straightforward, with 75 inserted in the appropriate 
places. 

Combining the n different Pqa with the possible index 
combinations of the Green's functions results in a total of 
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4, 12 and 12 (n = 1,2,3) choices to compute the vertex 
amplitude Ao(j>) in Eq. all of which should provide 
the same result for the final renormalisation constant in 
the absence of lattice artefacts. Because of the different 
sized momentum components in different lattice direc- 
tions, the expected discretisation errors vary depending 
on the directions selected by the indices of the projector. 
We reflect these artefacts coming from breaking contin- 
uum 0(4) symmetry to lattice hypercubic symmetry in 
the systematic error of our final results. With additional 
lattice spacings and the use of partially- twisted boundary 
conditions, we could eliminate hypercubic lattice arte- 
facts in the continuum limit 1 78], ual ■ 



3. Quark field renormalisation 

In general, the renormalisation condition Eq. (|5ip re- 
quires knowledge of the field renormalisation Z q to ob- 
tain Zq. However, in the present calculation only ratios 
of renormalisation factors of operators with one, two or 
no derivatives appear, Eqs. (|25p and (|2"8"1) . Combining 
this with our renormalisation condition leads to, 



{12} 
2 



{13} 
3 



{14} 
4 



{23} 
3 



{24} {34} 

4 4 



Zq 



n=2,3 



(m) 



a bare 



=i(p) 



(61) 



where the explicit Z q dependence drops out. As men- 
tioned earlier, we can use either the vector or axial- vector 
bilinears in this ratio thanks to chiral symmetry. We fol- 
low our earlier procedure [59[ and average A 7p and 
(Ay/A^ in the reference) to obtain our best answer. The 
analysis is also performed with A^ only and the differ- 
ence of the two enters our systematic error. 



4- Results for renormalisation factors 
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{12} {13} {14} {23} {24} {34} 
1112 2 3 

{123} {124} {124} {134} {234} {234} 
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0.672 
0.670 
0.668 
0.666 
0.664 
0.662 
0.660 
0.658 
0.656 
0.654 
0.652 



0.459 
0.458 
0.457 
0.456 
0.455 
0.454 
0.453 
0.452 
0.451 
0.450 

{123} {123} {124} {134} {134} {234} 
13 2 14 3 

FIG. 7. Results for A^ r n c =2 (A^ n e =3 ) on the top (bottom) for 
a fixed momentum (apf = 1.78201, p T = (2,2,3,8). The la- 
bels above and below the plots show the indices of the Green's 
function {pi} (top) and projector (bottom). The disagree- 
ment between the different projections is due to lattice arte- 
facts. 
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Compared to [59( the reduced statistical errors make 
previously hidden systematic effects apparent and quan- 
tifiable [63j and affect the way we extract the renormal- 
isation factors. We start by considering different pro- 
jectors for a fixed momentum p^ of the external quarks, 
see Fig. [71 The results should be independent of the ro- 
tation and size of the momentum components used for 
the projector. The smaller statistical errors now reveal 
a disagreement due to lattice artefacts. We combine all 
choices for our best answer and account for the spread 
in our systematic error, improving previous estimates. 
Our general recipe to obtain the renormalisation factors 
follows. The ratio of bare vertex amplitudes is extrap- 
olated linearly to the chiral limit m q —¥ —m res for each 
momentum. Only in the chiral limit can we remove the 
running of our data points and match them to a contin- 
uum scheme. So by using (6^, [Ml Iz3 we take our results 
from the RI'/MOM scheme at scale /i 2 = p 2 to a com- 
mon scale p? = 4 GeV 2 and convert to MS at that scale. 



The values thus obtained are then linearly interpolated to 
p 2 = (2 GeV) 2 within our momentum window to obtain 
Zo,n=2,3/Zo, n =i at a scale p = 2 GeV. 

The central value is computed from the averaged values 
from all projectors and index combinations. A standard 
bootstrap analysis provides the statistical error which is 
inflated with ^/x 2 /d.o.f. (the PDG scale- factor 80] ) from 
the interpolation. Several effects are taken into account 
for the systematic error. Lattice artefacts are the domi- 
nant effect. To estimate those, we perform the analysis 
for all projectors separately as indicated above and chose 
the highest and lowest result for each momentum for the 
interpolation. From the two fits, the larger deviation 
from the central value then constitutes the systematic 
error from discretisation effects (labelled 'spread' in the 
final table). This is a conservative approach for the dis- 
cretisation error. Taking random choices of projectors 
(or rather their direction) for each momentum and look- 
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ing at the la width of the range of results for many of 
those picks would lead to a smaller error. We account for 
missing higher order terms in the continuum perturbative 
calculation via the slope of the momentum interpolation, 
using the difference of our results at p 2 = (2 GeV) 2 and 
(OGeV) 2 , indicated by 'slope', as a measure. We note, 
however, that we cannot disentangle perturbative and 
discretisation errors here and thus double count some of 
the discretisation effects. Another source of systematic 
error is the strange quark mass, kept fixed at m s = 0.04 
in our simulation. We deal with that as described at the 
end of section IV. F in [59|, estimating an error from half 
the linear dependence (slope) multiplied by the strange 
quark mass, m s . This error is labelled 'Ara s '. The last 
contribution to the systematic error is from the chiral 
symmetry breaking evident when comparing our vector 
and axial- vector operators [5^, [73|, [8l[ . This is estimated 
by the difference of the final results when taking the axial- 
vector bilinear (n — 1) or the averaged vector and axial- 
vector bilinear for the ratio in Eq. ([Fiji (labelled 'V— j4'). 
Adding the four contributions in quadrature gives our 
systematic error. 

To illustrate some of the steps mentioned above, we 
include in Fig. [5] two examples of the extrapolation to 
the chiral limit. Shown are extrapolations for all our five 
momenta, for the renormalisation factors for one and two 
derivatives. In Fig. [3] we show the renormalisation factors 
before and after we remove the running, again for one 
and two derivatives. Once at a common scale, the data 
points arc much flatter indicating the validity of the scale 
conversion and the momentum window. Also included is 
the linear fit and our final result. 

Our final values for the ratios of renormalisation fac- 
tors are given in Table|Xj These have been obtained using 
vector like operators. Results from the axial- vector oper- 
ators are almost identical and show no low energy effects 
from breaking chiral symmetry, as for bilinears. We note 
that the renormalisation factors are significantly different 
from one and deviate substantially from the perturbative 
results. Thus nonperturbative renormalisation looks im- 
perative here. 



C. Renormalised Results 

We now use the renormalisation factors from the pre- 
vious section to convert our bare lattice results to MS 
at [i = 2 GeV. The local matrix elements in Eq. © re- 
quire the renormalisation factors Ze> {p(j} , Zdd.dd arid 
ZDD,dd as defined in Eqs. (1551) and (|55|) . The first two 
are computed nonperturbatively while for the last one we 
use the perturbative result. From Eq. (|3"5j) we see that 
the mixing term requires the computation of a matrix 
element with an operator insertion of Oqq. This is sim- 
plified since we use the ratios (|2"5|) and (|28p to extract 
the moments of the distribution amplitudes. Within the 
ratios, the matrix element with the operator Oqq differs 
from the denominator only by the momentum factors and 
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FIG. 8. Linear extrapolations of the renormalisation factors 
to the chiral limit. The top shows Zq 1 ^® m , the bottom 
plot is for 2 , Q I n / ^ J OM . The momenta are increasing from top 



to bottom and we have {ap) 2 
and 1.9801. 



1.2947, 1.4392, 1.6374, 1.7820 



thus does not have to be computed separately. It con- 
tributes a constant shift to the result. To summarise: 



MS 



Zq 



Z A 



^DD,dd 



Z A 



(62a) 
(62b) 



With our best nonperturbative results from Table [X] 
and the perturbative result for the mixing term from 
Table IIXI (computed at the same scale of fj, = 2 GeV, 
0.027(8)), we arrive at the renormalised mo- 



Z A 



ments of the distribution amplitudes given in Table IXII 



The contribution from the mixing term in Eq. (|62b|) is 
small so using the perturbative result for Zjjd^qq is not 
a drawback. Even if the correction of a nonperturba- 
tive result for ZoD,dd is as sizeable as for Zq {p ^ ] or 
Zdd,dd, the overall contribution remains comparable to 
our present error on Zjjd^dd/Za- Hence our results are 
essentially renormalised nonperturbatively. 
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TABLE X. Final results for the renormalisation factors in MS at p = 2 GeV. Results are given for both lattice sizes with all 
systematic errors. The perturbative results are also shown for comparison. 







Zos^x/Za 


Zdd, 


dd/Za 




16 3 x 32 


24 3 x 64 


16 3 x 32 


24 3 x 64 


central value 


1.54575 


1.52893 


2.06064 


2.02800 


statistical error 


0.00249 


0.00081 


0.00482 


0.00149 


spread 


0.02968 


0.01809 


0.03702 


0.01534 


slope 


0.00470 


0.00743 


0.00097 


0.02285 


Am s 


0.00089 


0.00232 


0.00469 


0.00992 


V -A 


0.00723 


0.00602 


0.00938 


0.00760 


total error 


0.03102 


0.02061 


0.03879 


0.03026 


best result 




1.5289(8) (206) 


2.028(1)(30) 


perturbative result 


1.24(3) 


1.45(5) 



TABLE XI. Final results in the chiral limit in MS at n = 2 GeV for both of our lattice volumes. Here the first error is statistical, 
the second includes systematic errors from m s , discretisation and renormalisation. 
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16 3 x 32 


0.25(1)(2) 


0.035(2) (2) 


0.25(1)(2) 


0.25(2)(2) 


0.037(1) (2) 


0.25(1)(2) 


0.24(1)(1) 


24 3 x 64 


0.28(F)(2) 


0.036(1)(2) 


0.26(1)(2) 


0.27(1) (2) 


0.043(2)(3) 


0.25(2)(2) 


0.25(2)(1) 



V. SUMMARY 

We have computed the first or first two lowest non- 
vanishing moments of the distribution amplitudes of the 
7T, K, K* , p and <f> mesons, using nonperturbative renor- 
malisation of the lattice operators, with final numbers 
given in Table IXD Apart from the uncertainty in m s for 
the first moments, systematic errors mainly come from 
the renormalisation procedure. Within the current sta- 
tistical errors on our data we do not see any finite size 
effects. With only one lattice spacing we can also only es- 
timate a formal discretisation error of 0(a 2 Aq CD ) « 4% 
from the 0(a)-improved DWF action and operators; this 
is included in our sytematic error. The result for {^)k 
in Table IXj s upercedes but is compatible with our earlier 
result in [4l|, [42j, which was obtained on the 16 3 x 32 
ensembles only and used perturbative renormalisation. 

Converting the lowest moment of the kaon distribu- 
tion amplitude to the first Gegenbauer moment ar K = 
0.061 (2) (4), we find it in agreement with sum rule results 
from Eq. (|16[) but with a much reduced uncertainty. We 
compare our results to those from the QCDSF Collabo- 
ration 23] in Table IXIII (preliminary results for the first 
moment of the vector meson distribution amplitudes are 
also available from QCDSF [27]). The results for (^} K 
differ significantly. However, we observe that our mea- 
surements correspond to pion masses in the range 330- 
670 MeV and are for 2 + 1 dynamical flavours, whereas 
the QCDSF results are for pion masses around 600 MeV 
and higher, with 2 dynamical flavours. For one data 
point from each collaboration where the pion and kaon 
masses are comparable, the (£, 1 )k values differ by about 
one standard deviation. These points occur for the small- 
est values of m 2 K — m 2 from each collaboration; for larger 



values the points, and therefore slopes in m 2 K — , differ. 
We plan to improve our results in the near future by re- 
ducing the systematic uncertainties. We will improve the 
nonperturbative calculation of the renormalisation fac- 
tors by including the total derivative mixing term. We 
will also have an additional lattice spacing allowing us to 
estimate the continuum results, including using partially- 
twisted boundary conditions to remove hypercubic lattice 
artefacts 0, |22|. Increased statistics on the 24 3 x 64 
lattice should also improve our conclusions about finite 
volume effects. 



ACKNOWLEDGMENTS 

The calculations reported here used the QCDOC com- 
puters [3 [82|, at Edinburgh University, Columbia 
University and Brookhaven National Laboratory (BNL). 
The Edinburgh QCDOC system was funded by PPARC 
JIF grant PPA/J/S/1998/00756 and operated through 
support from the Universities of Edinburgh, Southamp- 
ton and Wales Swansea, and from STFC grant 
PP/E006965/1. At BNL, the QCDOC computers of 
the RIKEN-BNL Research Center and the USQCD Col- 
laboration were used. The software used includes: 
Chroma H, QDP++ and the CPS QCD codes 13, 
supported in part by the USDOE SciDAC program; the 
BAGEL [52[ assembler kernel generator for many of the 
high-performance optimized kernels; and the UKHadron 
codes. We thank the University of Southampton for 
access to the Iridis computer system used in the cal- 
culations of the nonperturbative renormalisation factors 
(with support from STFC grant ST/H008888/1). DB, 
MAD, JMF, AJ, TDR and CTCS acknowledge support 



19 



TABLE XII. Comparison to other lattice results (both for MS at fi = 2GeV). 
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FIG. 9. These plots show the scale dependent Z's and Z's 
for a fixed scale [i — 2 GeV with the running successfully re- 
moved (both in RJ'/MOM). The top (bottom) plot is for the 
one (two) derivative case. Also included are the linear inter- 
polation to the final result with the statistical error indicated 
by the error band. Our result at fi = 2 GeV is then shown 
with error bars for the statistical and systematic errors. 
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